function z = fcov(X,d)
mu = mean(X,d);
X0 = bsxfun(@minus, X, mu);
z = X0*X0'/size(X,d);